import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
from math import log
from sklearn.feature_selection import SelectKBest, f_classif
from sklearn import preprocessing
from sklearn import metrics
import warnings
%matplotlib inline
import matplotlib.pyplot as plt
#defaults
plt.rcParams['figure.figsize'] = (20.0, 20.0)
plt.rcParams.update({'font.size': 10})
plt.rcParams['xtick.major.pad']='5'
plt.rcParams['ytick.major.pad']='5'
plt.style.use('ggplot')
warnings.filterwarnings("ignore", category=FutureWarning )
warnings.simplefilter(action='ignore', category=(UserWarning,RuntimeWarning,FutureWarning))
pd.set_option('display.max_columns', 500)
pd.set_option('display.max_rows', 500)
pd.options.mode.chained_assignment = None
For demonstration of open source Interpretable ML packages I have used UCI credit card default dataset.
UCI credit card default data: https://archive.ics.uci.edu/ml/datasets/default+of+credit+card+clients
The UCI credit card default data contains demographic and payment information about credit card customers in Taiwan in the year 2005. The data set contains 23 input variables:
LIMIT_BAL: Amount of given credit (NT dollar)
SEX: 1 = male; 2 = female
EDUCATION: 1 = graduate school; 2 = university; 3 = high school; 4 = others
MARRIAGE: 1 = married; 2 = single; 3 = others
AGE: Age in years
PAY_0, PAY_2 - PAY_6: History of past payment;
PAY_0 = the repayment status in September, 2005;
PAY_2 = the repayment status in August, 2005; ...; PAY_6 = the repayment status in April, 2005.
The measurement scale for the repayment status is: -1 = pay duly; 1 = payment delay for one month; 2 = payment delay for two months; ...; 8 = payment delay for eight months; 9 = payment delay for nine months and above.
BILL_AMT1 - BILL_AMT6: Amount of bill statement (NT dollar).
BILL_AMNT1 = amount of bill statement in September, 2005; BILL_AMT2 = amount of bill statement in August, 2005; ...; BILL_AMT6 = amount of bill statement in April, 2005.
PAY_AMT1 - PAY_AMT6: Amount of previous payment (NT dollar).
PAY_AMT1 = amount paid in September, 2005;
PAY_AMT2 = amount paid in August, 2005; ...; PAY_AMT6 = amount paid in April, 2005.
default payment next month: Target Variable - whether or not a customer defaulted on their credit card bill in late 2005. Default (1) and Non-Default (0)
Data Prepration and Feature Engineering¶
data = pd.read_excel('default of credit card clients.xls', skiprows=1)
data = data.rename(columns={'default payment next month': 'DEFAULT_NEXT_MONTH'})
data.head(2)
for var in ['PAY_0', 'PAY_2', 'PAY_3', 'PAY_4', 'PAY_5', 'PAY_6']:
data[var] = data[var].apply(lambda x: 0 if x<=0 else x)
# Caluclate average payment and log of payment
data['PAY_AMT_AVG'] = data[[col for col in data.columns if col.startswith('PAY_AMT') ]].mean(axis=1)
data['PAY_AMT_AVG_log'] = data['PAY_AMT_AVG'].apply(lambda x: log(x+1))
#
for i in np.arange(1,7):
data['PAY_REL_AMT_'+str(i)] = data['PAY_AMT'+str(i)]/data['PAY_AMT_AVG']
#
# log of payments
for i in np.arange(1,7):
data['PAY_REL_AMT_log_'+str(i)] = data['PAY_AMT'+str(i)].apply(lambda x: log(x+1))
#Feature Engineering of bill amount
data['BILL_AMT_AVG'] = data[[col for col in data.columns if col.startswith('BILL_AMT')]].mean(axis=1)
data['BILL_AMT_AVG_log'] = data['BILL_AMT_AVG'].apply(lambda x: log(x+1) if x>0 else 0)
#
# bill sign as a separate feature
for i in np.arange(1,7):
data['BILL_AMT_SIGN_'+str(i)] = data['BILL_AMT'+str(i)].apply(lambda x: float(x>0))
#
#log of credit limit
data['LIMIT_BAL_log'] = data['LIMIT_BAL'].apply(lambda x: log(x+1))
data['LIMIT_BAL_CAT'] = pd.cut(data['LIMIT_BAL'], range(0, np.max(data['LIMIT_BAL']), 10000), right=False)
#
data['SEX'] = data['SEX'].astype('category').cat.rename_categories(['M', 'F'])
data['MARRIAGE'] = data['MARRIAGE'].astype('category').cat.rename_categories(['NA', 'MARRIED', 'SINGLE', 'OTHER'])
data['AGE_CAT'] = pd.cut(data['AGE'], range(0, 100, 10), right=False)
education_dict = {0:'other', 1:'graduate school', 2:'university', 3:'high school', 4:'other', 5:'other', 6:'other'}
data['EDUCATION'] = data['EDUCATION'].apply(lambda i: education_dict[i])
data = pd.get_dummies(data, columns=[col for col in data.columns if data[col].dtypes not in ['int64', 'float64']])
data.head(2)
#Imputing the missing values
for name in data.columns:
if pd.isnull(data[name]).sum() > 0:
median = data[name].median()
data[name] = data[name].apply(lambda x: median if pd.isnull(x) else x)
#
y = data['DEFAULT_NEXT_MONTH']
X = data.drop(['DEFAULT_NEXT_MONTH', 'ID'], axis = 1)
#
selector = SelectKBest(f_classif, 30)
selector.fit(X, y)
top_indices = np.nan_to_num(selector.scores_).argsort()[-30:][::-1]
selector.scores_[top_indices]
X_prep = X[X.columns[top_indices]]
scaler = preprocessing.MinMaxScaler()
X_transformed = scaler.fit(X_prep).transform(X_prep)
X_transformed_df = pd.DataFrame(X_transformed,columns=X.columns[top_indices])
X_train, X_test, y_train, y_test = train_test_split(X_transformed_df, y, test_size=0.2, random_state=1)
X_train.head(2)
Open Source Techniques for Interpretable Machine Learning¶
1. Permutation Feature Importance (PFI)¶
(Global Feature Importance Method)¶
PFI is an algorithm that computes importance scores for each of the feature variables of the dataset. The importance measures are determined by computing the senstivity of a model to random permutations of feature values.
In other words an importance score quantifies the contribution of a certain feature to the predictive performance of a model in terms of how much a choosen evaluation metric deviates after permuting the values of that feature.
Intuition behind a permutation importance is that a feature is “important” if altering or permuting its values increases the model error, because in this case the model relied on the feature for the prediction. A feature is “unimportant” if altering or permuting or shuffling its values leaves the model error unchanged, because in this case the model ignored the feature for the prediction.
Permutation Feature Importance provides global insight into model behavior.
from IPython.display import Image
Image("PFI.png")
# Train Random Forest Classifier
from sklearn.ensemble import RandomForestClassifier
rf_model = RandomForestClassifier().fit(X_train, y_train)
y_pred_rf = rf_model.predict_proba(X_test)
#
fpr, tpr, thresholds = metrics.roc_curve(y_test, y_pred_rf[:,1:2])
print('RF AUC on Test Data', metrics.auc(fpr, tpr))
# Using ELI5 Calculate Feature Importance using Permutation Importance
import eli5
from eli5.sklearn import PermutationImportance
perm = PermutationImportance(rf_model,scoring = 'roc_auc').fit(X_test, y_test)
eli5.show_weights(perm, top=100, feature_names = X_test.columns.tolist())
Interpreting Permutation Importance¶
The values towards the top are the most important features, and those towards the bottom matter least. The first number in each row shows how much model performance decreased with a random shuffling. The number after the ± measures how performance varied from one-reshuffling to the next.
Some time we see negative values for permutation importances. In those cases, the predictions on the shuffled (or noisy) data happened to be more accurate than the real data. This happens when the feature didn't matter (should have had an importance close to 0), but random chance caused the predictions on shuffled data to be more accurate.
2. LIME - Local Interpretable Model-Agnostic Explanations¶
LIME (Ribeiro et. al. 2016) is an implmentation of local surrogate model and is used to explain individual predictions of black box models. Surrogate models are trained to approximate the predictions of the underlying black box model.
The idea is quite intuitive. First, forget about the training data and imagine you only have the black box model where you can input data points and get the predictions of the model. You can probe the box as often as you want. Your goal is to understand why the machine learning model made a certain prediction. LIME tests what happens to the predictions when you give variations of your data into the machine learning model. LIME generates a new dataset consisting of permuted samples and the corresponding predictions of the black box model. On this new dataset LIME then trains an interpretable model, which is weighted by the proximity of the sampled instances to the instance of interest. The interpretable model can be anything for example Lasso or a decision tree. The learned model should be a good approximation of the machine learning model predictions locally, but it does not have to be a good global approximation. This kind of accuracy is also called local fidelity.
Algorithm:
Generate a fake dataset from the example we’re going to explain.
Use black-box estimator to get target values for each example in a generated dataset (e.g. class probabilities).
Train a new white-box estimator, using generated dataset and generated labels as training data. It means we’re trying to create an estimator which works the same as a black-box estimator, but which is easier to inspect. It doesn’t have to work well globally, but it must approximate the black-box model well in the area close to the original example.
To express “area close to the original example” user must provide a distance/similarity metric for examples in a generated dataset. Then training data is weighted according to a distance from the original example - the further is example, the less it affects weights of a white-box estimator.
Explain the original example through weights of this white-box estimator instead.
Prediction quality of a white-box classifer shows how well it approximates the black-box classifier. If the quality is low then explanation shouldn’t be trusted.
Example of LIME technique using LIME Package:
from xgboost import XGBClassifier
xgb_model = XGBClassifier()
xgb_model.fit(X_train.values, y_train)
#
y_pred_xgb = xgb_model.predict_proba(X_test.values)
#
fpr, tpr, thresholds = metrics.roc_curve(y_test, y_pred_xgb[:,1:2])
print('XGB AUC on Test Data', metrics.auc(fpr, tpr))
#all categorical columns
categorical_cols = [col for col in X_test.columns if X_test[col].dtypes not in ['int64', 'float64']]
#all feature names
all_cols = X_test.columns
import lime
from lime.lime_tabular import LimeTabularExplainer
explainer = lime.lime_tabular.LimeTabularExplainer(X_train.values ,class_names=[0, 1], feature_names = all_cols,
categorical_names=categorical_cols, mode='classification')
y_test[:5]
y_pred_xgb[:,1:2][:5] #probabilities for default(1) are shown
predict_fn = lambda x: xgb_model.predict_proba(x)
exp = explainer.explain_instance(X_test.values[3], predict_fn, num_features=10)
exp.show_in_notebook()
predict_fn = lambda x: xgb_model.predict_proba(x)
exp = explainer.explain_instance(X_test.values[4], predict_fn, num_features=10)
exp.show_in_notebook()
3. SHAP¶
SHAP – SHapley Additive exPlanations – explains the output of any machine learning model using Shapley values. Shapley values have been introduced in game theory since 1953 but only recently they have been used in the feature importance context. SHAP belongs to the family of “additive feature attribution methods”. This means that SHAP assigns a value to each feature for each prediction, the higher the value, the larger the feature’s attribution to the specific prediction. It also means that the sum of these values should be close to the original model prediction.
SHAP unified several existing feature attribution methods such as LIME, Deep Explainer and more and it theoretically guarantees that SHAP is the only additive feature attribution method with three desirable properties:
Local accuracy: The explanations are truthfully explaining the ML model
Missingness: Missing features have no attributed impact to the model predictions
Consistency: Consistency with human intuition (more technically, consistency states that if a model changes so that some input’s contribution increases or stays the same regardless of the other inputs, that input’s attribution should not decrease)
import shap
(A) SHAP for Tree Based Models (Xgboost)¶
from xgboost import XGBClassifier
xgb_model = XGBClassifier()
xgb_model.fit(X_train.values, y_train)
#
y_pred_xgb = xgb_model.predict_proba(X_test.values)
#
fpr, tpr, thresholds = metrics.roc_curve(y_test, y_pred_xgb[:,1:2])
print('XGB AUC on Test Data', metrics.auc(fpr, tpr))
print(y_pred_xgb[:,1:2][:5]) #First 5 predicted rows
print(y_test[:5]) #First 5 actual values
For 4th Observation Actual Outcome = 1, Predicted Probability = 0.74137545¶
# load JS visualization code to notebook
shap.initjs()
explainer = shap.TreeExplainer(xgb_model)
shap_values = explainer.shap_values(X_test.values)
# visualize the first prediction's explanation
shap.force_plot(explainer.expected_value, shap_values[3,:], X_test.iloc[3,:])
# SHAP shows log of odds in XGBOOST, so converting back to probability
1/(1+np.exp(-(explainer.expected_value + sum(shap_values[3,:]))))
explainer.expected_value, shap_values[3,:], X_test.iloc[3,:]
For 5th Observation Actual Outcome = 0, Predicted Probability = 0.34667933¶
# visualize the first prediction's explanation
shap.force_plot(explainer.expected_value, shap_values[4,:], X_test.iloc[4,:])
# SHAP shows log of odds in XGBOOST, so converting back to probability
1/(1+np.exp(-(explainer.expected_value + sum(shap_values[4,:]))))
explainer.expected_value, shap_values[4,:], X_test.iloc[4,:]
Global Feature Importance using SHAP¶
To get an overview of which features are most important for a model we can plot the SHAP values of every feature. The plot below sorts features by the sum of SHAP value magnitudes over all samples, and uses SHAP values to show the distribution of the impacts each feature has on the model output. The color represents the feature value (red high, blue low).
# summarize the effects of all the features
shap_values = explainer.shap_values(X_train.values)
shap.summary_plot(shap_values, X_train)
Mean absolute value of the SHAP values for each feature
shap.summary_plot(shap_values, X_train, plot_type="bar")
(B) SHAP with DeepExplainer for Neural Network based models (Keras and tensorflow)¶
from keras.models import Sequential
from keras.layers.core import Dense, Activation, Dropout
input_dim = X_train.shape[1]
nn_model = Sequential()
nn_model.add(Dense(30, input_shape=(input_dim,), activation='relu'))
nn_model.add(Dense(15, activation='relu'))
nn_model.add(Dense(1, activation='sigmoid'))
nn_model.compile(optimizer='adam', loss='binary_crossentropy', metrics=['accuracy'])
#
history = nn_model.fit(np.array(X_train), np.array(y_train),
batch_size=25, epochs=5, verbose=2,
validation_split=0.2)
score = nn_model.evaluate(np.array(X_test), np.array(y_test), verbose=0)
print('Test log loss:', score[0])
print('Test accuracy:', score[1])
nn_y_pred = nn_model.predict_on_batch(np.array(X_test))[:,0]
fpr, tpr, thresholds = metrics.roc_curve(y_test, nn_y_pred)
print('Neural Network AUC', metrics.auc(fpr, tpr))
print(nn_y_pred[:5]) #First 5 predicted values
print(y_test[:5]) #First 5 actual values
# select a set of background examples to take an expectation over
background = X_train.values[np.random.choice(X_train.shape[0], 100, replace=False)]
# explain predictions of the model
e = shap.DeepExplainer(nn_model, background)
For 4th Observation Actual Outcome = 1, Predicted Probability = 0.6936611¶
# load JS visualization code to notebook
shap.initjs()
shap.force_plot(e.expected_value, e.shap_values(X_test.values[3].reshape(1,-1))[0], X_test.iloc[3,:])
e.expected_value, e.shap_values(X_test.values[3].reshape(1,-1))[0], X_test.iloc[3,:]
For 5th Observation Actual Outcome = 0, Predicted Probability = 0.4703175¶
# load JS visualization code to notebook
shap.initjs()
shap.force_plot(e.expected_value, e.shap_values(X_test.values[4].reshape(1,-1))[0], X_test.iloc[4,:])
e.expected_value, e.shap_values(X_test.values[4].reshape(1,-1))[0], X_test.iloc[4,:]
(C) SHAP Model agnostic with KernelExplainer (explains any function)¶
from sklearn import linear_model
# Create logistic regression object
model_logistic = linear_model.LogisticRegression()
# Train the model using the training sets
model_logistic.fit(X_train, y_train)
#
y_pred_regr = model_logistic.predict_proba(X_test.values)
#
fpr, tpr, thresholds = metrics.roc_curve(y_test, y_pred_regr[:,1:2])
print('Logistic Regression AUC on Test Data', metrics.auc(fpr, tpr))
# use Kernel SHAP to explain test set predictions
X_train_summary = shap.kmeans(X_train, 10)
explainer = shap.KernelExplainer(model_logistic.predict_proba, X_train_summary)
# As per shap documentation if training size is big we have to cluster using kmeans
# rather than use the whole training set to estimate expected values, we summarize with
shap_values = explainer.shap_values(X_test, nsamples=100)
# plot the SHAP values for the Setosa output of the first instance
shap.force_plot(explainer.expected_value[0], shap_values[0][0,:], X_test.iloc[0,:])
explainer.expected_value[1]
#Actual values first five
y_test[:5]
#Predicted Probabilites by logsitic regression of 3rd, 4th, and 5th observations
y_pred_regr[:5][:,1:2]
For 4th Observation Actual Outcome = 1, Predicted Probability = 0.76505556¶
# plot the SHAP values
shap.force_plot(explainer.expected_value[1], shap_values[1][3,:], X_test.iloc[3,:])
explainer.expected_value[1] + sum(shap_values[1][3,:])
explainer.expected_value[1], shap_values[1][3,:], X_test.iloc[3,:]
For 5th Observation Actual Outcome = 0, Predicted Probability = 0.35220043¶
# plot the SHAP values
shap.force_plot(explainer.expected_value[1], shap_values[1][4,:], X_test.iloc[4,:])
explainer.expected_value[1], shap_values[1][4,:], X_test.iloc[4,:]
explainer.expected_value[1] + sum(shap_values[1][4,:])
Probabilities by logistic regression - coff and interept
model_logistic.coef_
model_logistic.intercept_
import math
def sigmoid(x):
return 1 / (1 + math.exp(-x))
#probabilites from coeff and intercept
[sigmoid(x) for x in (np.dot(X_test[:5].values, model_logistic.coef_.T ) + model_logistic.intercept_).reshape(1,-1)[0]]
4. LOCO - Leave-One-Covariate-Out¶
LOCO creates local interpretations for each row in a training or unlabeled score set by scoring the row of data once and then again for each input variable (e.g., covariate) in the row. In each additional scoring run, one input variable is set to missing, zero, its mean value, or another appropriate value for leaving it out of the prediction. The input variable with the largest absolute impact on the prediction for that row is taken to be the most important variable for that row’s prediction. Variables can also be ranked by their impact on the prediction on a per-row basis.
Its the most simplest of all and can deteriorate in accuracy when complex nonlinear dependencies exist in a model.
LOCO is the most simplest and seems to be not a good technique and can be implemented in straight forward way as shown below:
X_test_5_df = X_test[:5].copy()
#Original Predition for top 5 rows in test data
y_pred_5 = xgb_model.predict_proba(X_test_5_df.values)[:,1:2]
y_pred_5
LOCO_df = pd.DataFrame(columns = X_test_5_df.columns)
# Lets suppose if we want to get the feature attribution for top 5 observations in test data
for i in (X_test_5_df.columns):
# train and predict with Xi set to missing
X_test_5_df_cpy = X_test_5_df.copy()
X_test_5_df_cpy[i] = np.nan
#print(X_test_5_df_cpy.head())
preds_y_test_5_df = xgb_model.predict_proba(X_test_5_df_cpy.values)[:,1:2]
#print(preds_y_test_5_df)
print('LOCO Progress: ' + i)
# subtract the LOCO prediction from original prediction
LOCO_df[i] = (y_pred_5 - preds_y_test_5_df).reshape(1,-1)[0]
LOCO_df
X_test_5_df
References:¶
Permutation Feature Importance:
https://christophm.github.io/interpretable-ml-book/feature-importance.html
https://eli5.readthedocs.io/en/latest/blackbox/permutation_importance.html
LIME:
https://eli5.readthedocs.io/en/latest/blackbox/lime.html
https://christophm.github.io/interpretable-ml-book/lime.html
https://github.com/marcotcr/lime
SHAP
https://github.com/slundberg/shap
LOCO
https://github.com/h2oai/mli-resources/blob/master/notebooks/loco.ipynb